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APPLICATION OF CONTINUED FRACTIONS FOR FAST EVALUATION 
OF CERTAIN FUNCTIONS ON A DIGITAL COMPUTER 


Amnon Bracha 


Abstract 

The purpose of this paper is to develop a method for evaluation 
of certain elementary functions on a digital computer by the use of 
continued fractions. The time required for this evaluation is drastically 
reduced by using "short" operations like shift and add, instead of 
multiplications. Functional consistency is the most important factor that 
allows the expansion of a function into a continued fraction. Several cases 
are discussed and in particular the solution of the quadratic equation 


is discussed in more detail to demonstrate the convergence of the method. 


1. INTRODUCTION 

The idea of using continued fraction representations for generating 
a solution to a limited class of quadratics was first introduced by 
Robertson [4]. 


Consider the finite continued fraction with k partial numerators 


“i 
ei” 


D; and k partial denominators ee 1s Bn soon Is WaOSS “Welle is i526 5 
(a 
Me = PL 
foe ch ya 
=a 
as, 


A convenient way of writing (1.1) is 


ere ee Pk. | 
aren pat hs ! 


AL. and B. are determined from the recursions: 


Be Cg ea | 
(1.2) 


with initial values: 


> 
| 
© 
> 
I] 
i 
fH 


O al 


Bo = db BL 


OG 
iG aif} @lleeno wihehy AL and B. can be separately and simultaneously determined 
in two binary arithmetic units in k-1 addition times if the Ps and qd, are 
chosen to be simple in the binary sense. 

The digit set for Py and q. for the purposes of this paper is 


(U2 Lhe it will be stated later that the continued fraction ** 


Be 


assumes all values in the limit, over the interval 


The range defined in (1.3) includes the range of normalized floating 
point binary fractional parts, [2 /2,A J, This property indicates that a suitabli 
continued fraction representation exists, such that conversion to conventional — 
binary can be achieved by repetitive use of two binary adders in parallel, 


followed by a division to determine the quotient a 
k 


The main reason for selecting Pi» a.e{1, 1/2), i=1, 2, »-. is\ theta 


four multiplicative operations required for each iteration in (1.2) are 


3) 
reduced to "shift" and "add" operations. These operations will be called 
"short" operations throughout this paper, mainly because the time required 
to perform these operations is shorter than the time required to perform 
"long" operations, e.g., Multiplication, Division. 

The purpose of this paper is to develop algorithms for fast 
evaluation of certain elementary functions by using "short" operations in 
several registers simultaneously. In order to be able to do so we make 


use of functional consistency which will be defined at the end of section 2. 


Determination of selection rules for p and gq in each iteration 
is an important step for the development of the algorithm. Selection 
rules were extensively studied by Trivedi [5], where a complete set of 
such rules were developed for the quadratic equation. The set of 
selection rules that is used in this paper is described in section 3. 

In section } we generalize our results to a higher degree 
polynomial and in section 5 we show two more cases where our analysis 


is applicable. 


2. BILINEAR TRANSFORMATIONS AND THE RICCATI EQUATION. 


in this seetion we develop a special case or tierandeisy cmc 
of Wynn [8]. 
The general continued fraction will be regarded as a 
sequence of bilinear transformations of the form: 
Px 


(Ao) f= ——=— , ke=l1,2,... 
k ati ay ? 


where f, (x) is a function of x. Therefore 


Pp Pp p 
Ge) t= = =... —— 
om ED Ga ntl 
A +f A 
ia faded jad 
(2.3) = ane eB ee Ie Jalint=turcl eeh es r P aphe hs 


a oie 


where the functions A and Bo Satisfy the recursion 


a A +p A 
n cial n-1 Pa n-2 
(2,1) 
B =| = 
Tae ae) Se ee A 
with the initial values 
Kg © AL TS 
© le aaa 
Bo = B = 
ate ee coil 


For the purposes of this paper Ps ds» ty 2 aay ehdlll lee siuLacimacl 
from the digit set {1/2, 1} so that the recursion (2.4) can be performed 


by using only "short" operations. 


WN 


The main purpose of this section is to show that there exist 
funetions for which bilinear transformations of the form (2.1) can be 
used, and such that the functions te k=i1, 2, .«.-, are consistent, 


1.6. Equal. 


Consider the Riccati equation 


2 
(2.5) ee ne he See. © 


where y, is a function of the variable x, and ay, by and c, are Punetlons! Ot: 
x or constants. The property of this equation as noted by P. Wynn [7] is 
that if the dependent variable vy is replaced by the bilinear 

transformation (2.1) then the functions te iS dig 25 ooo BlUSO Senoilsiive eae 


Riccati equation 
2 
u _ 
(2.6) eee by, 2c = 0 


We develop below the recursion for the coefficients of the (Eanes equation 


by means of the coefficients of the ue equation. 


2 
' _ 
Let v7 ay, tf bY, oe Cc. = '© 
benune Keth Riceati equation. 
(2.1) = / 
From (2.1) we have y= So eps te. (ll. 2) 
k q+ Tey es Waele 
-P,, Vy, 
then since ie =e 
(ie) 
we have 
1h. Nise P P 
etal i k es k Bee 


If we multiply by (96%) to normalize the coefficient of Tea we get 
k . 


2 

e 2cmd: Cag 
! k 2 k-k k°k 
USL” i. “eal (b+ ) Vagy 7 (OP, HO, at 


y= 6 
Py Py 


and the recursion that follows is 


ie eee 
ar db PL 
2E, G 
k*k 
(2.70) Bag = Pe = 5 
k 
ec 
k*k 
Seg, OT SP oh. = a 


We note that all the operations involved in (2.7) require only "short" 


CHAomeAialeas , Sines jochela Py and q,. are Simple binary constants. 


Lemma 1: [1] 


2] 
A = bd = ha, Cy = constant ke leer oe 


Proof: We use the recursion (2.7) and get 


22 
of fre eee HOC Oh : anor oe 
lend PV er ied mele p OES Siri 
k 19 
k 
2 ® 
Py neta 
. Py > 2 Bh bes sc ots 


We define below the term functional consistency. 


Definition: For a substitution of the form 


(x) “ 

0S) see CE 1, 2, «ses 

k ke] x 

where f(x) is a function of x, p, and q, are constants; if f(x) and 
Fie (&) is the same function, then we have functional consistency for 


(x), 


3, SOLUTION OF ax~ + bx - c =0 


We show now how the solution of a quadratic equation 
with two distinct roots of opposite sign, and in particular 
the square root problem, can be found by the technique of 
section 2. 


Let 


Ged) Ae bp = © =O 


be a given quadratic equation. 


The substitution we use is of the form 


(3.2) =! 
202 Gr 
Gla 
where Ps, a € (1, 1/2) i 2 il, 2, 
Hor. thie Re step we have 
po Pp 
k k 
a Re % ee ae = @ = © k = Ween 
oh ret 
or 
io a (Ae e =1, 1.) ede: i= pe = 18,19, 6 0 
eager Tee ey ee, Ieee “ete 7 SRI 
The recursion that follows is 
aren) ie 
(5-5) Big 7 Ee eee 


k =a Oe 


Q 
os 
+ 
fe 

H] 

i 

Q 

Pa 
a 

+ 

xe) 
be 

+ 

io 
Ts 
a 

Ke) 
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and the resulting quadratic equation is 


} 
(3.4) Peet cn * Mey een ~ Sea = ° 


This method of approximating the solution of (3.1) can 


be used if we develop a technique for selecting Py and qd. Kira 25 


from the coefficients of the k-th quadratic equation, i.e., ay» by 
and Che 

The following lemma will be stated without a proof. 
Lemma 2: 


Let the k@? approximation to the solution of (4.1) be 


M= max xy = Jo 
Ti <= Mulia oe = va-1 
iL 2 


Using the result of the Lemma it can be seen that the functional consistency 


of the procedure can be achieved in each step if 
(3:5) meh Sw ote Sal Oty ae 


By imposing condition (3.5) we need only one set of 


selection rules for Py and qd hea 2... LOmmrner ance: jm iM): 


We develop now a set of selection rules for P,, and Che 
Io tS, 2 cong woe We cwracimenle equesiem [5]. 
We write below a version of (3.1). 


Let 


(5.6) 4 ESS 


where it is assumed that c. 7 0, b. 20. (4 


iL L 5 EO, eialel jm Se = iN 


a 
We will find Py and dy such that 


(Bn) m= 


where mS x, SMandp,, q, € fi, ale). 


Clearly, we have four possibilities, and for cach patron 
CliLae} : 
Py and qd, we feta ditferent Xp 
We take now the inverse approach. We assume that condition 


(4.5) exists for x., and find the range for x, for each pair of p 
2 1 


tL 


and q,. We start with a pair p, =1 andq, = 1/2. From (3.7) we have: 


For x,, =e (lower bound for x, ) 


1 e(eNe-1) _ 
(3.8) xX, 2m SE CF 0522, 
1 i/e+e ii 
and for i> (Vo-1)/2 (upper bound for x, ) 
(3.9) a hoe SN 2 ih, 


The result is that for x, in the range defined by (3.8)-(3.9), 


al 


we can choose Py = 1 and qq = 1/2. Since x) is the unknown we use 


(4.6) in order to find the allowable range for p, =l anda, = WD: 


We have 


> 1 > 
(3.10) Tay) ea eae mere =) 06922): 


oy ea 


Since (3.10) is possible for any x, in the range (3.8)-(3.9), we 


A 
conclude that for the range 


«0.508, € @ & Ne es 


(ea) 0. 522b 1 8 = fl Hi 


A 


we choose Py = i and: q, = nH / 2 
Similarly we write below the ranges for each of the remaining 
Possibilities: 


For 


G12) (J2-1)b, + We-1)°a, en a(ve-1)b, +4W2-1)*a,, choose 


2 
avo-1 eNo-1 NE i 
(3.13) cau Dam er aceae) a, Soy <b, + 5a,, choose 


ale 


(3.14) ——b_+ ( 


2 
es (J2-1)b, + W2-1) a,» choose 


Ww 
i) 
fe 
I 


Mfr 


oe > 4, =1. 


The result is that the entire range is divided ante four 
sections, (4.11)-(3.14), and for each section we can choose a pair 
of p, and q, and such that condition (4.5) for X, be satisfied. 
Clearly, if we have to do two multiplications for each 
selection range in order to find Py and qj» our procedure is) anerivemendce 
In the analysis that follows we make use of an important feature of the 
ranges defined in (3.11)-(3.14); this is the existence of overlapping 
between any two consecutive ranges. 
This means that in the overlap regions we have a freedom 
of selecting the pairs Py and qa between two sets of such constants. 
We will use this freedom in order to simplify our selection algorithm 
by defining a line of selection inside the overlap region and such 
that the coefficients of by and ay Will -be simple in the binary sense: 
Before we define the selection lines, we note that rave yor 
convergence of the method was found to be strongly dependent on these 


lines. The first set of selection Lines can be the upper ines mam 


each range, for 


aL 
Ga) 0.4 Ub, < 0.1713a, then p, = 57% Ns 
(525) a = OGOyb, = 0.5 then Ugo: vee) Queer 
; ith : lama vee SL iL Dyed aed Bye 
oS 0. 828b,, < 0. 686a,, then p, =1,4 =1; 


otherwise Py = dks qd, = = : 


Experimentally this set of rules gave the best rate of convergence. 


To simplify the constants that appear in (3.15) we use 


simple binary constants with at most two non-zero binary digits, 


We have, 
Ge Ons7ob. = O.t5625e Ghene | Da x am all 
ale E fla aS 1 all Cie, Selle Z 
O7e25p. = 0 then oe 2, 
mer oh = Os 98, lle 2a Kae a 
(4.16) is ae 
e, - 0.75b, = 0.6254, then py =2 50, =1; 


otherwise Pi = Je qd, = 5 5 


These selection rules involve only short operations. ‘The procedure 
can be carried out now in the same way for Xn and so on. 

The algorithm described in this section involves the following 
steps, for the k iterations: 


(1) Equation (4.4) is given, then use selection 


rules (3.16) to find p, and q,. 


(2) Use the results of step (1) and iterate on 
(aoe 


(4) Use the recursion (4.4) and find (3.4) for k+1. 


(4) Check if A,/B,. reached the required precision. 
This check can be done only once if the number 
of jterations required to achieve certain 
precision is known. The analysis below for the 
rate of convergence gives the necessary information 
to find such numbers. If the required precision 


is not reached proceed to step (1) for k+l. 


Two examples are shown in appendix A. The rules of selections 


which were used are (3.16). 


1 


The following theorem assures convergence of the algorithm. 


Theorem lL: 


Let 


| A 2 n 
od De Me Saree eral Women a 
(3-17) ile B 1 ate do aun Gee qa ap xe n=) (2 eerie 


be the exact solution of (3.1), which was found by substitution (4.2), 
with the set of sellection rules (3.16) and recursion (4.3) for the 


coefficients of (4.1). 
A 
n th ‘ : 
Let Bp? be the n~ approximation to x 
n 
e-Q, arise N sien tae stow alll im = N, 


1 Then for every 


Ricoik ek NPR, OOSEGIEWE TwhayeG Acie x, SAGLSLLSS wae ee equation Gale then 


x, is a solution of (S61). 


We will study now a relation between om and § because of 


=D” 


a well known property of odd and even convergents of a continued fraction. 


Define 


Mee 0) aka Sreme ets 


then we only have to prove that Th <i? 


We have [2] (Chapter 1): 


nt+e2 
A/B - AB GL Sao Be eat! © Ba 
a IVS aN B 7% ra 
fall “ah 2 (21) PyPor**Py_ (Int In Pe /B Be 
Fray tel : By-2 


Get Sa Pat a 


Since Pal r Sait. a P,/X, *: In? In+1 = 1 and B, F Ber mel y Pee A 


we conclude that: 


7 PA *n nee 
i B Zi. B 

D +a, n=J 1+ Tab 

oe n Ee 


All the quantities in Th are positive, Le Sil, charck ‘slate 
result follows. 
We will use now the analysis of Theorem 1 to study the rate 


of convergence, by finding an upper bound to rn 


tn 


1 -—x 


Max T = Max 
n 
1+ 


be 8 


Therefore, 


ia che ee Nes 
B ie eae Te 
n-2 
Finally we have, 
5 _ Bel 5 
a ee es 
Max TS ma Phledd eNES 0.2929 , 
1 v2 eee 
= “ 


Therefore for n sufficiently large, the error oe is reduced by 
a factor which is less than or equal to 0.2929 for each pair of 


additional iterations. 


4. SOLUTION OF A HIGHER DEGREE POLYNOMIAL 


We show now how one solution of the cubic equation 
can be found by the method of section 43. 


Let 


2 e oe 
(4.1) an; + bx, + cx, - d, =0 
be a given cubic equation. We use the substitution (3.2) and we 


get for the yon step: 


He Pr Py 
* 'b =o) 
is Cla (eee) oes. * 


or 


oe (Baie. C, Pye) 44 % (Bd).ay- 2c, Pye ,.-b,, Pi) 41 


2 
(4,07 +b, PE +0, Dd - d,.a° ) ma 


The recursion relations between the coefficients of the 


ie weunie equation and the Goa cubic equation are therefore: 


Teel 
qd > th Pi 
SS ee = Behe Ge-2C4P cP, Pi 
qa =a,,Det Pedy oP, 4-4, a7 Is tly, ome ec 


the resulting cubie equation is 


3 2 : 
Sea iea | Perse) Seen Se 


18 


For the selection rules we use an analysis similar to 
that of section 43. First observe that the bounds given in 
(3.8)-(3.9) for the case p = 1 and q =1/2 are valid. Therefore 
we can write an expression, similar to (3.10) for the cubic 


equation: 


IV 


(25) 2 — > _ = 0.520 


The result is that for 


0, 5220, + 0.520°b. + 0.522%a. < d s 2c, AB) & avea, 


dL dL JL iL 


we choose p, =1 anda, = /2> 


For the remaining cases we have 


(V2-1)e,+W2-1)*,+ (/2-1)7a, sd, s2W2-1)c,+h W2-1)%b, +8 (/2-1)7a, 


? 


choose Pa = JL 5 qd, = 4s 


ovo-1 eve-1\* eVo-1\ ? Jo 1 Jo 
ec, + |———] b a <d, <s—eca, +=b, + — 
ion 7 1 ji > Sa pr a ae ee Sama 
, choose es WB, qd, = Ly Bs 
Jo-1 NiBai\\ Jo-1\? ke 3 
a Sh + 3 Das a, 2 = (Ve-1)c, & @2-1) By * (Vo-1) as. 


choose P) = L/D. = a 


hh 


We are ready to write now a set of selection rules similar 


HO) (Go ALS) « i.€., the upper dine in each range will bevour seileciion 


line: 


- ) - — — ha 

d, - O.4lke, - 0.1713b, $ 0.071a, then p, 1/2, q, = 4: 

(4.4) d, - 0.707c, - 0.5b, = 0.35358, then p, = w/e), q, = 1/2; 
d, - 0,828c, - 0.686b, < 0. 5683a, then p, =1, 4, =413 


otherwise B= as qd = ie : 


We note that as in (3.16) the constants which appear in 
(4.4) can be simplified. 

For the proof of convergence and the rate of convergence 
we can use the analysis of section 3, and therefore we develop a 
method to approximate one positive solution of a cubic equation. 

The procedure can be generalized now to higher degree 
polynomials. The necessary steps are 

(1) to write the recursion for the coefficients 

of the polynomial 


(2) to develop the selection rules by using an 
argument similar to (3.10) and (4.3) 


(3) to simplify the coefficients in the selection 


rules. 


The result is an algorithm which always converges to one 
positive solution. 
In appendix B we give two examples for approximation 


of a positive root of a cubic equation. 


20 


5. CONTINUED FRACTION EXPANSIONS OF CERTAIN FUNCTIONS BY THER USE 
OF THE GENERAL SOLUTION OF THE RICCATI EQUATION 
Let 
; 2 
(5) AL) y' + ay’ + by + c =O 


be a Riccati equation, with a, b, and c constants. 


In order to find the solution of (5.1) we integrate by 


parts: 
di 
(G2) : aii 
: ay +by+e 
and the solutions are: 
dy iL eaytb- Vbo-hac 
(5) oR Te eee LS eee eee 
ay +by+e J p°-hac eay+bt vp2-hac 
= = Arctgh a when Be ohae > Ox 
nb -hac Jpo-hac 
=2 e ‘ 
= biDay 9 when b -lac = O53 
e2aytb 
= ss arctg — » when po -hac SO), 
Vhac-b Vac-b 


The above solutions can be used now for the continued 
fraction expansion of the inverse functions which appear explicitly 
in the solution. 


We start with the case 


(Ses) Vi = we xx 


The Riccati equation for (5.4) is 


We note that -A =Vhac-b2 =2. 

Now we use a bilinear transformation of the form (2.1). 
The result is a differential equation of the type (2.6) with the 
recursion (2.7). By Lemma 1 it follows that the solution for each 
equation k is of the type (5.3) with-A> 0, and therefore we get for 
the «oh step: 


2a, td, 


-x,-, STarece 5 
where d,. is a constant of integration. 


The solution is therefore 


Except for the first part of the solution which is a linear 
transformation, we see the consistency of the method, because if 

a set of selection rules are developed for tg x it can be used for 
each step and therefore evaluation of this function will be possible. 


Another important function which can be included is oe 


We have 
ya eye 
The ku step solution is 
2a, y, tb, -1 

Page gaa 2a, yj, tb Fl ? 

or 
b tl 1 
Mice fou er eae 


Again we note that if a set of selection rules can 
be developed for e* then it is possible to carry the process Sere 
each step and therefore to find the continued fraction expansion 
for the exponential function. 


For the case where A = 0 we have several possibilities: 


(a) b=a=O, y' + ¢ =O with the 


SO ULURvLCia, We = Hee te El 
2 
(o) bse =0, y' + ay =O with the 
solution 
eel 
Vu ~ Oxtb 


(o) 0, 1 SO, @ 4 0, sad i alee oeO 
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